#include <stdio.h>
#include <stdlib.h>
#include "linsyst.h"
#include "logtime.h"

#define LOGFILE "mjr-cuda.log"

#define NBLOCK 64
#define NUMTHREAD 512

void transformaestrela(syst * lsys);
void jacobi_richardson(syst * lsys);

int main(int argc, char** argv){
	int i, j;
	float teste, el;
	syst * lsys;

	if(argc != 2){
		printf("Usage: %s FILE\n\n", argv[0]);
		return 1;
	}

	initlogtime();

	// Carrega arquivo
	lsys = loadsyst(argv[1]);
	if(lsys == NULL) return 1;

	regtime(LOADFILE);

	//Aplica o metodo
	transformaestrela(lsys);
	jacobi_richardson(lsys);

	setit(lsys->max_it);


	//Imprime resultado
	printf("Iterations: %d\n", lsys->max_it);
	
	i = lsys->ref;
	el = lsys->el_digref;
	teste = 0;
	for(j = 0; j < lsys->size; j++){
		if(j != i) teste += MAT_A(lsys,i,j)*(lsys->vetX0[j])*el;
		else teste += (lsys->vetX0[j])*el;
	}

	printf("RowTest: %d => [%f] =? %f\n", i, teste, lsys->vetB[i]*el);
	destroysyst(lsys);


	//Registra execucao no arquivo de log
	savelog(argv[1], LOGFILE);

	return 0;
}


void transformaestrela(syst * lsys){
	int i, j;
	float el_dig;

	for( i = 0; i < lsys->size; i++){
		el_dig = MAT_A(lsys,i,i);
		MAT_A(lsys,i,i) = 0;
		for( j = 0; j < lsys->size; j++)
			MAT_A(lsys,i,j) /= el_dig;

		lsys->vetB[i] /= el_dig;
	}
}


/* reducao que soma um vetor em cuda*/
__device__
void sum_tree_like_reduction_sum(float * vector, int vector_size, float *sum){
	int j, k;

	for(k = vector_size/2; k > 0; k >>= 1){
		__syncthreads();
		for( j = threadIdx.x; j < k; j += blockDim.x)
			vector[j] += vector[k+j];
	}

	__syncthreads();
	if(threadIdx.x == 0) *sum = vector[0];
}


/* calcula novo X*/
__global__
void calc_Xk(int size, float * mA, float * vB, float * vX0, float * vX1)
{
	__shared__ float accum_threads[NUMTHREAD];
	int first, last, i, j, k;
	float A_X0_parc_product;

	for(i = blockIdx.x; i < size; i += gridDim.x)
	{
		first = i*size; //inicio da linha
		last = first + size; //fim da lina
		A_X0_parc_product = 0.0;
		k = threadIdx.x; //numero da thread no bloco

		//realiza a soma parcial da thread
		for( j = k + first; j < last; j += blockDim.x, k += blockDim.x)
			A_X0_parc_product += mA[j]*vX0[k];

		accum_threads[threadIdx.x] = A_X0_parc_product; //copia para o vetor shared do block
		
		__syncthreads();
		sum_tree_like_reduction_sum(accum_threads, NUMTHREAD, &A_X0_parc_product); //reducao soma

		if(threadIdx.x == 0)
			vX1[i] = vB[i] - A_X0_parc_product; // Salva elemento de X calculado pelo bloco
	}
}

/*Sum tree modificada para registrar o maximo*/
__device__
void sum_tree_like_reduction_max(float * vector, int vector_size, float *sum){
	int j, k;

	for(k = vector_size/2; k > 0; k >>= 1){
		__syncthreads();
		for( j = threadIdx.x; j < k; j += blockDim.x)
			if(vector[k+j] > vector[j]) vector[j] = vector[k+j];
	}

	__syncthreads();
	if(threadIdx.x == 0) *sum = vector[0];
}

// variavel no device para registrar o erro
__device__ float xerro = 0.0;

//Copia X1 para X0 e calcula o erro da iteracao
__global__
void update_X0(int size, float * vX0, float * vX1)
{
	__shared__ float max_x_per_th[NUMTHREAD];
	__shared__ float max_dif_per_th[NUMTHREAD];
	int i;
	float dif, max_x, max_dif;

	max_x = 0.0; //maior x[i] encontrado pela thread
	max_dif = 0.0; //maior diferenca encontrada pela thread

	for(i = threadIdx.x; i < size; i += blockDim.x){
		//atualiza maior x[i]
		if( vX1[i] > max_x) max_x = vX1[i];
		else if( -vX1[i] > max_x) max_x = -vX1[i];

		//atualiza maior diferenca
		dif = vX1[i] - vX0[i];

		if(dif > max_dif) max_dif = dif;
		else if(-dif > max_dif) max_dif = -dif;

		// capia X1[i] para X0[i]
		vX0[i] = vX1[i];
	}

	//Copia maximos parciais no vetor compartilhado
	max_x_per_th[threadIdx.x] = max_x;
	max_dif_per_th[threadIdx.x] = max_dif;
	__syncthreads();

	sum_tree_like_reduction_max(max_x_per_th, NUMTHREAD, &max_x); //reducao maximo de x
	sum_tree_like_reduction_max(max_dif_per_th, NUMTHREAD, &max_dif); //reducao maximo da diferenca

	if(threadIdx.x == 0 && max_x > 0.0)	xerro = max_dif/max_x; //calcula erro
}

void jacobi_richardson(syst * lsys){
	int it;
	int size;
	float *devA, *devB, *devX0, *devX1;
	float herro;

	size = lsys->size;

	//Aloca memoria do dispositivo
	cudaMalloc( (void**)&devA, size * size * sizeof(float));
    cudaMalloc( (void**)&devB, size * sizeof(float));
    cudaMalloc( (void**)&devX0, size * sizeof(float));
	cudaMalloc( (void**)&devX1, size * sizeof(float));


	//Copia A* B* para memoria do dispositivo
	cudaMemcpy( devA, lsys->matA, size * size * sizeof(float), cudaMemcpyHostToDevice);
	cudaMemcpy( devB, lsys->vetB, size * sizeof(float), cudaMemcpyHostToDevice);
	cudaMemcpy( devX0, lsys->vetB, size * sizeof(float), cudaMemcpyHostToDevice);

	//registra tempo para transferencia
	regtime(MCOPY_HD);

	//executa iteracoes
	for(it = 0; it < lsys->max_it; it++){

		//calcula novo x
		calc_Xk<<<NBLOCK,NUMTHREAD>>>(size, devA, devB, devX0, devX1);
		update_X0<<<1,NUMTHREAD>>>(size, devX0, devX1);

		//atualiza x antico e calcula erro
		cudaMemcpyFromSymbol(&herro, "xerro", sizeof(xerro), 0, cudaMemcpyDeviceToHost);
		if(herro < lsys->prec){
			it++;
			break;
		}
	}

	//registra tempo de execucao do metodo
	regtime(METJACOBI);


	//copia resultado para memoria do host
	cudaMemcpy( lsys->vetX0, devX0, size * sizeof(float), cudaMemcpyDeviceToHost);

	//registra tempo para copiar	
	regtime(MCOPY_DH);

	//libera memoria do dispositivo
	cudaFree( devA);
    cudaFree( devB);
    cudaFree( devX0);
	cudaFree( devX1);

	//registra numero de iteracoes
	lsys->max_it = it;
}



